home *** CD-ROM | disk | FTP | other *** search
/ Shareware Grab Bag / Shareware Grab Bag.iso / 007 / c_math.arc / JN.C < prev    next >
Encoding:
C/C++ Source or Header  |  1983-07-02  |  1.8 KB  |  108 lines
  1. /*
  2.     floating point Bessel's function of
  3.     the first and second kinds and of
  4.     integer order.
  5.  
  6.     int n;
  7.     double x;
  8.     jn(n,x);
  9.  
  10.     returns the value of Jn(x) for all
  11.     integer values of n and all real values
  12.     of x.
  13.  
  14.     There are no error returns.
  15.     Calls j0, j1.
  16.  
  17.     For n=0, j0(x) is called,
  18.     for n=1, j1(x) is called,
  19.     for n<x, forward recursion us used starting
  20.     from values of j0(x) and j1(x).
  21.     for n>x, a continued fraction approximation to
  22.     j(n,x)/j(n-1,x) is evaluated and then backward
  23.     recursion is used starting from a supposed value
  24.     for j(n,x). The resulting value of j(0,x) is
  25.     compared with the actual value to correct the
  26.     supposed value of j(n,x).
  27.  
  28.     yn(n,x) is similar in all respects, except
  29.     that forward recursion is used for all
  30.     values of n>1.
  31. */
  32.  
  33. #include <math.h>
  34. #include <errno.h>
  35.  
  36. int    errno;
  37.  
  38. double
  39. jn(n,x) int n; double x;{
  40.     int i;
  41.     double a, b, temp;
  42.     double xsq, t;
  43.     double j0(), j1();
  44.  
  45.     if(n<0){
  46.         n = -n;
  47.         x = -x;
  48.     }
  49.     if(n==0) return(j0(x));
  50.     if(n==1) return(j1(x));
  51.     if(x == 0.) return(0.);
  52.     if(n>x) goto recurs;
  53.  
  54.     a = j0(x);
  55.     b = j1(x);
  56.     for(i=1;i<n;i++){
  57.         temp = b;
  58.         b = (2.*i/x)*b - a;
  59.         a = temp;
  60.     }
  61.     return(b);
  62.  
  63. recurs:
  64.     xsq = x*x;
  65.     for(t=0,i=n+16;i>n;i--){
  66.         t = xsq/(2.*i - t);
  67.     }
  68.     t = x/(2.*n-t);
  69.  
  70.     a = t;
  71.     b = 1;
  72.     for(i=n-1;i>0;i--){
  73.         temp = b;
  74.         b = (2.*i/x)*b - a;
  75.         a = temp;
  76.     }
  77.     return(t*j0(x)/b);
  78. }
  79.  
  80. double
  81. yn(n,x) int n; double x;{
  82.     int i;
  83.     int sign;
  84.     double a, b, temp;
  85.     double y0(), y1();
  86.  
  87.     if (x <= 0) {
  88.         errno = EDOM;
  89.         return(-HUGE);
  90.     }
  91.     sign = 1;
  92.     if(n<0){
  93.         n = -n;
  94.         if(n%2 == 1) sign = -1;
  95.     }
  96.     if(n==0) return(y0(x));
  97.     if(n==1) return(sign*y1(x));
  98.  
  99.     a = y0(x);
  100.     b = y1(x);
  101.     for(i=1;i<n;i++){
  102.         temp = b;
  103.         b = (2.*i/x)*b - a;
  104.         a = temp;
  105.     }
  106.     return(sign*b);
  107. }
  108.